home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 007 / fourier.bas < prev    next >
BASIC Source File  |  1985-04-15  |  4KB  |  124 lines
  1. 2 '***********************************************************************
  2. 4 '*     FOURIER SMOOTHING WITHOUT THE FAST FOURIER TRANSFORM PROGRAM    *
  3. 6 '*             By Eric E. Aubanel and Keith B. Oldham                  *
  4. 8 '***********************************************************************
  5. 10 CLS
  6. 12 INPUT "ENTER NUMBER OF DATA POINTS";N
  7. 14 REM LEAVING R AND I ARRAYS UNDIMENSIONED LIMITS VALID VALUES OF E TO <=10
  8. 16 DIM X(N),X1(N),U(N),V(N)
  9. 18 FOR I=0 TO N-1
  10. 20      INPUT "ENTER DATAPOINT VALUE";X(I)
  11. 22 NEXT I
  12. 24 GOSUB 60
  13. 26 PRINT "THE SMOOTHED DATA VALUES ARE:"
  14. 28 FOR I=0 TO N-1
  15. 30      PRINT "X(";I+1;") = ";X1(I)
  16. 32 NEXT I
  17. 34 END
  18. 36 REM FOURIER ALGORITHM SUBROUTINE BEGINS AT LINE 60.  LINE NUMBERS ARE THE           SAME AS FOR THE HP VERSION OF THE SUBROUTINE
  19. 60 PI=3.141593
  20. 70 PRINT "NUMBER OF TRANSFORM POINTS TO BE KEPT";
  21. 80 INPUT E
  22. 90 IF E>INT((N+1)/2) THEN PRINT "E TOO LARGE":GOTO 70
  23. 100 IF E<>INT(E) OR E<=1 THEN GOTO 70
  24. 110 IF E<=Q THEN 870
  25. 120 REM
  26. 130 IF Q<>0 THEN 330
  27. 240 'CALCULATE R(0)
  28. 250 G=0
  29. 260 FOR J=0 TO N-1
  30. 280 G=G+X(J)
  31. 290 NEXT J
  32. 300 R(0)=G/N
  33. 310 Q=1
  34. 320 REM
  35. 330 PRINT "WORKING ON R(K) TRANSFORM CALCULATIONS"
  36. 340 J2=INT((N-1)/2)
  37. 350 P1=INT(LOG(2*J2-1)/LOG(2))
  38. 360 FOR K=Q TO E-1
  39. 370     J1=J2
  40. 380     S=PI*K*2/N
  41. 390     C=COS(S):S=SIN(S)
  42. 400     FOR J=1 TO J1
  43. 410             L=2*J-1
  44. 420             U(J)=X(L)*C+X(L+1)
  45. 430             V(J)=X(L)*S
  46. 440     NEXT J
  47. 450     S=2*S*C:C=2*C*C-1
  48. 460     FOR P=1 TO P1
  49. 470             U(J1+1)=0:V(J1+1)=0
  50. 480             J1=INT((J1+1)/2)
  51. 490             FOR J=1 TO J1
  52. 500                     L=2*J-1
  53. 510                     U=U(L)*C-V(L)*S+U(L+1)
  54. 520                     V(J)=U(L)*S+V(L)*C+V(L+1)
  55. 530                     U(J)=U
  56. 540             NEXT J
  57. 550             S=2*S*C:C=2*C*C-1
  58. 560     NEXT P
  59. 570     R(K)=(X(0)+(U(1)*C+V(1)*S))/N
  60. 580 NEXT K
  61. 590 REM
  62. 600 PRINT "WORKING ON I(K) TRANSFORM CALCULATIONS"
  63. 610 FOR K=Q TO E-1
  64. 620     J1=J2
  65. 630     S=2*PI*K/N
  66. 640     C=COS(S):S=SIN(S)
  67. 650     FOR J=1 TO J1
  68. 660             L=2*J-1
  69. 670             U(J)=-(X(L)*S)
  70. 680             V(J)=X(L)*C+X(L+1)
  71. 690     NEXT J
  72. 700     S=2*S*C:C=2*C*C-1
  73. 710     FOR P=1 TO P1
  74. 720             U(J1+1)=0:V(J1+1)=0
  75. 730             J1=INT((J1+1)/2)
  76. 740             FOR J=1 TO J1
  77. 750                     L=2*J-1
  78. 760                     U=U(L)*C-V(L)*S+U(L+1)
  79. 770                     V(J)=U(L)*S+V(L)*C+V(L+1)
  80. 780                     U(J)=U
  81. 790             NEXT J
  82. 800             S=2*S*C:C=2*C*C-1
  83. 810     NEXT P
  84. 820     I(K)=-((U(1)*C+V(1)*S)/N)
  85. 830 NEXT K
  86. 840 REM
  87. 850 IF E>Q THEN Q=E
  88. 860 REM
  89. 870 PRINT "WORKING ON INVERSE TRANSFORM"
  90. 880 REM
  91. 890 'CALCULATE X1(0)
  92. 900 F1=0:F2=0
  93. 910 FOR K=1 TO E-1
  94. 920     T=R(K)
  95. 930     F1=F1+T
  96. 940     F2=F2+K*K*T
  97. 950 NEXT K
  98. 960 X1(0)=R(0)+2*(F1-F2*(1/E/E))
  99. 980 REM
  100. 990 P1=INT(LOG(2*E-3)/LOG(2))
  101. 1000 FOR J=1 TO N-1
  102. 1010    T2=E*E
  103. 1020    FOR K=1 TO E-1
  104. 1030            F=1-K*K/T2
  105. 1040            U(K)=R(K)*F:V(K)=-(I(K)*F)
  106. 1050    NEXT K
  107. 1060    K1=E-1
  108. 1070    S=2*PI*J/N
  109. 1080    C=COS(S):S=SIN(S)
  110. 1090    FOR P=1 TO P1
  111. 1100            U(K1+1)=0:V(K1+1)=0
  112. 1110            K1=INT((K1+1)/2)
  113. 1120            FOR K=1 TO K1
  114. 1130                    L=2*K-1
  115. 1140                    U=U(L)*C-V(L)*S+U(L+1)
  116. 1150                    V(K)=U(L)*S+V(L)*C+V(L+1)
  117. 1160                    U(K)=U
  118. 1170            NEXT K
  119. 1180            S=2*S*C:C=2*C*C-1
  120. 1190    NEXT P
  121. 1200    X1(J)=R(0)+2*(U(1)*C+V(1)*S)
  122. 1220 NEXT J
  123. 1230 RETURN
  124.